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Abstract  The  solution  of  linear  systems  having  circulant  coefficient  matrices  is  considered  in  this 
paper.  This  kind  of  systems  occur  in  many  applications:  prediction,  time  series  analysis,  spline 
approximation,  difference  solution  of  partial  differential  equations,  etc.  The  methods  presented 
here  are  more  efficient  than  the  Toeplitz  type  methods  and  are  ba?ed  on  the  fast  Fourier  transform 
as  well  as  the  circulant  factorization  of  the  ^banded  circulant”™'-*" — - 
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1.  Introduction 

Many  problems  in  mathematics  and  applied  science  lead  to  the  solution  of  linear  systems  having 
circulant  coefficient  matrices,  which  are  related  to  the  periodicity  of  the  problems.  Examples 
are  prediction  [17],  time  series  analysis[l2],  spline  approximation[l],  solution  of  certain  partial 
differential  equations[l8],  and  possibly  in  many  other  applications. 

Circulant  matrices  belong  to  the  class  of  Toeplitz  matrices.  Linear  systems  with  Toeplitz 
matrices  can  be  solved  with  0(n3)  arithmetic  operations  [2,  11,  15,  19].  If  we  use  any  of  the 
methods  for  solving  Toeplitz  systems  to  solve  linear  systems  having  circulant  coefficient  matrices 
we  will  not  be  able  to  make  full  use  of  the  properties  of  the  circulant  matrices,  which  are  a)  the 
orthonormal  eigenvectors  of  all  n-by-n  circulant  matrices  are  the  columns  of  the  Fourier  matrix, 
the  elements  of  which  can  be  expressed  with  exponential  or  trigonometric  functions.  ?nd  6)  their 
eigenvalues  are  the  polynomials  of  the  primitive  n-th  root  of  unity  with  the  elements  of  the  matrices 
as  its  coefficients[9].  Such  nice  properties  make  it  simple  to  estimate  the  invertibility  of  the  matrix, 
to  compute  their  inverses(if  they  sxist) ,  and  to  solve  circulant  linear  systems.  Furthermore,  the  fast 
Fourier  transform  may  be  employed  to  calculate  the  eigenvalues  and  to  solve  the  systems  since  the 
components  of  the  eigenvectors  are  exponential  functions.  The  methods  presented  in  this  paper  are 
based  on  the  use  of  fast  Fourier  transform.  The  dominant  work  of  the  algorithms  is  in  performing 
the  fast  Fourier  transform,  and  therefore  the  algorithms  solve  circulant  systems  in  0(nlog2n) 
operations  as  opposed  to  the  0(n 2)  arithmetic  operations  required  by  using  other  algorithms  for 
solving  linear  systems  having  Toeplitz  coefficient  matrices. 

In  some  applications,  for  example,  in  the  finite  difference  solution  of  one  dimensional  elliptic 
equations  subject  to  periodic  boundary  conditions[5,  18]  and  approximation  of  periodic  functions 
using  splines[l],  “banded  circulant”  matrices  are  encountered.  Under  certain  coditions  such  matri¬ 
ces  can  be  factored  as  a  product  of  tow  simpler  circulants,  and  the  systems  may  then  be  solved  by 
using  the  Sherman-Morrison  formula,  or  its  block  version,  the  Woodbury  formula[l3].  This  method 
is  quite  competitive  with  Gaussian  elimination  both  in  terms  of  arithmetic  operations  and  storage 
requirements. 

In  the  case  where  multi-dimensional  problems  are  concerned  the  matrices  of  coefficients  of  the 
resulting  linear  systems  are  block  circulant  matrices.  After  some  transformations  and  permutations 
we  are  led  to  a  block  diagonal  matrix  with  circulant  blocks  on  the  diagonal.  This  reduces  the 
problem  to  the  solution  of  n  circulant  linear  systems,  which  may  be  performed  in  parallel.  An 
important  example  is  the  finite  difference  approximate  solution  of  elliptic  equation  over  a  rectangle 
with  periodic  boundary  conditions  [5,  18]. 
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In  §2,  we  develop  the  Fourier-circulant  method  for  solving  circulant  linear  systems.  In  §3  the 
method  for  decomposing  banded  circulant  matrices  into  two  circulants  and  solving  the  system  will 
be  described.  The  methods  presented  in  §2  and  §3  will  be  extended  to  block  circulant  systems  in  §4. 
§5  applies  the  results  of  §2,  §3  and  §4  to  elliptic  equations,  and  finally,  some  numerical  experiments 
will  be  given  in  §6. 

In  what  follows,  we  employ  the  notation  diag(6\ ,  <52, . . . ,  6n)  to  mean  the  n-by-n  diagonal  matrix 
with  diagonal  elements  6i ,  62, . . . ,  6n,  A*  the  complex  conjugate  transpose  of  the  matrix  A,  and  AT 
the  transpose  of  matrix  A. 


2.  Fourier- Circulant  Method 

Consider  the  system  of  linear  equations 


Cx  =  b. 


where  C  is  a  circulant  matrix  of  order  n,  and  x  and  b  are  n-vectors.  An  n-by-n  matrix  C  =  (c,y)  is 
a  circulant  matrix  if  c,y  =  Ci+ij+i,  and  the  subscripts  are  taken  modulo  n.  Thus  a  circulant  matrix 
can  be  written 


co 

Cl 

C2 

.  . .  C„_l 

Cn-l 

Co 

Cl 

• • •  Cn—2 

cn—2 

Cn— 1 

CO 

•  •  •  Cn_3 

| 

• 

; 

: 

Cl 

C2 

C3 

...  Co 

It  is  clear  that  C  contains  at  most  n  distinct  elements,  and  therefore  often  denoted  by 
(2.2)  C  —  circ(co,  ci,  C2,  •  •  • ,  cn— i). 

It  is  well  known[9]  that  if  C  is  a  circulant,  then 

(2-3)  C  =  Fdia<7(A1,A2,..„An)F*, 


where  F,  called  Fourier  matrix,  is  n-by-n  unitary  matrix  of  the  eigenvectors  of  matrix  C  with 
(i,  j)-element  n'suf''1^'1)  ,  and 


At  =  k  =  1, 2, . . . ,  n 


are  the  eigenvalues  of  C,  and  w  is  the  primitive  n-th  root  of  unity. 
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If  C  is  nonsingular,  then  the  inverse  of  C  is  also  a  circulant,  and  it  is,  from  (2.3),  given  by 


(2.5)  C'1  =  F  diag( Af1,  A2 1 , . . . , A ^)F'. 

Thus  the  solution  to  (2.1)  is 

(2.6)  x=F diag{\il , Aj  \ . . . , X~1)F*b. 

If  C  is  singular,  in  which  case  some  of  the  eigenvalues  of  C  are  zero,  the  inverse  of  C  will  not 
exist,  but  there  will  exist  an  unique  matrix  called  the  Moore-Penrose  generalized  inverse  of  C  and 
denoted  by  C+,  which  is  geven  by 

(2.7)  C+  =  Pdiaff(A+,A+,...,A+)F*, 


where 

(2.8) 


if  A*  =  0, 
if  At  ^  0, 


k=  l,2,...,n, 


and  the  least  squares  solution  of  (2.1)  is  given  by 


(2.9) 


x  =  C+b , 


which  is  the  solution  to  (2.1)  with  smallest  2-norm  when  the  system  is  consistent[3,  14]. 

It  is  easily  seen  that  premultiplying  a  vector  by  matrices  F*  and  F  may  be  accomplished 
by  the  fast  Fourier  transform(FFT)  and  its  inverse,  respectively,  and  the  eigenvalues  of  C  can  be 
calculated  via  FFT  [8].  Thus  the  algorithm  proceeds  as  follows. 

Algorithm  CIRS  (CIRculant  Solver)  solves  the  circulant  system  (2.1). 

1.  Transform  b  =  F*b  by  FFT. 

2.  Compute  the  eigenvalues  of  C  via  FFT. 

3.  Calculate  5  =  diag( Aj",A2  ,...,A+)6. 

4.  Tramsform  the  vector  5  to  obtain  the  solution  vector  x  =  F  h. 

endalgorithm 

The  algorithm  uses  the  fast  Fourier  tansform  three  times  and  is,  therefore,  an  O(nlog2n) 
algorithm. 


3.  Banded  Circulant  Systems 


In  the  case  where  only  cq ,  cj  and  c„_i  are  nonzero,  the  matrix,  denoted  by  Ct,  is  circulant 
tridiagonal ,  and  its  eigenvalues  can  easily  be  computed.  Assume  that  Ct  is  real,  then  the  eigenvalues 
are  given  by 


27T  27 r 

(3.1)  A*  =  Co  +  (ci  +  c„_i]  cos(A:  -  1) - h  i[ci  -  c„_i]  sin(A:  -  1)  — ,  k  =  1,2 , . . . ,  n, 

b  n 


where  i  =  \/— T.  If  the  matrix  is  real  and  symmetric(denoted  by  C4<),  i.e.  ci  =  cn_i,  then  the 
eigenvalues  are  real  and  given  by 

,  2;r 


(3.2) 


A*  =  co  +  2ci  cos(A:  —  1)- 


*  =  1,2,. 


In  these  situations,  the  method  presented  in  the  previous  section  may  be  employed.  Since  the 
matrices  Ct  and  Cst  have  the  cyclic  tridiagonal  structure,  we  are  able  to  take  advantage  of  the 
fact  that  under  certain  conditions  this  special  kind  of  matrices  can  be  factored  into  the  product 
of  a  circulant  lower  bidiagonal  matrix  and  a  circulant  upper  bidiagonal  matrix,  and  the  systems 
then  may  be  solved  by  using  the  Sherman-Morrison  formula[l3|.  In  this  section  we  will  develop  an 
algorithm  for  solving  the  systems 


(3.3) 

and 


Ctx  =  6, 


(3.4) 


C,tx  =  b, 


and  then  extend  it  to  general  banded  circulant  systems. 

We  assume,  without  loss  of  generality,  that  co  >  0.  It  was  shown  [6]  that  if  the  matrix  Ct  is 
strictly  diagonal  dominant,  then  there  exist  real  numbers  a,  0  and  7,  such  that  |/?|  <  1,  j-fj  <  1, 
and 


(3.5) 

where 


Ct  =  aLU , 


0  1 


L  = 


0\ 


0  1 
0  1 J 


It  is  easy  to  verify  that 


1  ' 

*  =2  [■ 


/*=^, 

a 

a 


1  1  \ 

1  1 

1  1 

1  1/ 
[co  +  Jcq-  4cicn-i  , 


To  solve  (3.3)  we  first  solve 

(3.7) 
and  then 

(3.8) 


Lh  =  a  16, 


Ux=h . 


where  h  is  an  auxiliary  vector  and  its  components  are  denoted  by  hj,  j  =  1, 2, . . . ,  n. 

Writing 

L  —  L  +  /JejeJ 

and 

0  =  U  +  -ye„e[, 

and  using  the  Sherman-Morrison  formula,  (3.7)  and  (3.8)  are  reduced  to  the  following  bidiagonal 


systems 

(3.9) 

Lh  =  a-l\b-  {fi ,(-/?)' 

{  i=i 

(3.10) 

Ux  =  h  —  ri^hii- 7),_1  e, 

.  1=1 

where 

'  1 

0  1 


0  1 

0  1. 


Cj  and  e„  are  the  first  and  the  last  column  of  the  n-th  order  identity  matrix,  respectively,  and 


(3.11) 

?  =  0/[l  -  (-/?)"], 

(3.12) 

»?  =  tt/[i  —  ( — y)"]. 

which  may  be  solved  by  forward  and  backward  sweeps,  respectively,  and  the  procedure  is  stable 
since  |/?|  <  1  and  |-y |  <  1.  The  algorithm,  thus,  may  be  formulated  as  follows. 

Algorithm  CTS(Circulant  Tridiagonal  Solver)  solves  circulant  tridiagonal  system  (3.3). 

1.  Compute  the  quantities  a,  0,  and  7  via  (3.6). 

2.  Compute  f  and  rj. 

3.  Solve  equation  (3.9)  and  (3.10). 

endalgorithm 

The  algorithm  solves  equation  (3.3)  with  0(5n)  arithmetic  operations,  which  is  more  favorable 
than  those  proposed  by  Bjorck  and  Golub  [4]  requiring  0(9n)  operations. 

The  method  studied  above  may  be  extended  to  “ banded  circulanf  matrices.  Consider  the 
symmetric  “banded  circulant”  of  the  form 


TO 


Assume  that  p  <§C  n,  the  order  of  the  matrix  C b.  If  Cb  is  positive  definite  and  elliptic[lO],  then 

there  exists  a  real  polynomial  l(z)  —  0o  + - h  0Pzp,  0o  >  0,  with  no  root  inside  the  unit  circle, 

such  that  the  characteristic  function  <f>(z)  of  A ,  defined  by 


4>{z)  =  cpzp  +  ■  ■  ■  +  cq  +  •  •  •  +  cpz  p, 


can  be  factored  as  4>{z)  =  l(z)  ■ 1(1/ z ).  There  are  several  methods  to  compute  this  factor(see  [10]). 
We  now  assume  that  the  factor  l(z)  has  been  computed.  It  is  easy  to  verify  that 


(3.13) 

where 


c6  =  llt, 


/  Po 


L  = 


0p 


0i  \ 


0p 


0p 


0o  J 


Cbx  =  d 


The  system 

(3.14) 

can  be  solved  by  solving  following  two  systems 

(3.15)  Ly  =  d 
and 

(3.16)  lTx  =  y. 


Let  E\  and  Ep  denote  the  n-by-p  matrices  consisting  of  the  first  p  columns  and  the  last  p 
columns  of  the  n-th  order  identity  matrix,  respectively,  and 

(0o 


L  = 


0p 


0p 


00  J 


Then  L  can  be  written 


(3.17) 


L  =  L  +  EiRE*, 


where 


(  ftp 


R  = 


fti\ 


V 


ftp) 


and  using  the  Woodbury  formula  the  inverse  is  given  by 


L~x  =  L~x  -  L~1Ei(R~1  +  El L~x Ei)~x Ej L~l, 


and  thus 


(3.18) 


y  =  L~ld  -  L~xEi(R~x  +  E% L~x Ei)~l E* IT1  d 
—  h  —  Wg, 


where  h,  W  and  g  are  the  solution  of  the  following  equations, respectively 


(3.19) 


Lh  =  d, 


(3.20) 


LW  =  Ei, 


(3.21) 


Bg  =  z, 


and 


2  —  (hrj-p+l  i  •  •  •  5  h n)T 

is  the  p-vector  with  the  last  p  elements  of  the  vector  h  as  its  components,  and 


(3.22) 


B  =  R~l  +  E^L-'Ei. 


To  compute  the  p-th  order  matrix  B,  we  first  solve  the  equation  (3.20).  Since  W  is  the  first  p 
columns  of  L~l,  and  L  is  a  lower  triangular  Toeplitz  matrix  and  so  is  its  inverse,  IV  is  defined  by 
the  first  column  of  L~l,  which  is  the  solution  of  the  equation 


(3.23) 


Lw=  (1,0,...  ,0)r, 
8 


and  can  be  computed  with  0(np)  operations. 

Denote  by  w\ ,  W2, . . . ,  wn  the  components  of  vector  w,  then  we  have 


(3.24) 


U>2 


•  Wn  Wn—  i  • 


(3.25) 


EjW  = 


Wn— p+l  Wn—p 


Wn-2p+2 


W„  Wn-i  ■  ■  ■  Wn-p+ 1 


which  is  a  Toeplitz  matrix  too. 

The  matrix  f?-1  is  an  upper  triangular  Toeplitz  matrix  and  can  be  calculated  with  0(p 2) 
operations.  Thus  B  is  Toeplitz,  so  solving  (3.21)  will  cost  0(p2)  operations.  Having  computed  B 
and  solved  the  equations  (3.19),  (3.20)  and  (3.21),  the  auxiliary  vector  y  can  be  found,  and  then 
we  can  solve  equation  (3.16)  in  a  similar  way.  Since 


l~T  =  L~t  -  L-TEpB~TETxL-\ 


the  solution  vector  x  is  given  by 
(3.26) 


x  =  r  —  Vs, 


where  r  =  (ri,  rj, . . . ,  rn)T  is  the  solution  of  the  equation 


(3.27) 


LT  r  =  y, 


t^n-p+l  p+ 2  Wn 

Wn—p  Wn—p+i  Wn—i 


(3.28) 


Wn-p+ i 


•  *'•  W  ^  ->  W"-  .*  .**.-.'*  .’ll-.  >  -  -  -  •  w  .  *  .  *  .  *  -  *  .  V  *.*  V 


Sr 


and  s  is  the  solution  of  the  equation 


(3.29)  BTs=(rur2,...,rp)T. 

The  asymptotic  operation  counts  of  the  method  would  be  0(5pn)  excluding  the  amount  of  work  to 
calculate  the  factor  l(z).  In  most  usual  case,  p  =  1  or  2,  and  finding  l(z)  does  not  cost  much  work. 
The  algorithm  may  be  summarized  as  follows. 

Algorithm  BCS  (Banded  Circulant  Solver)  solves  banded  circulant  system  (3.14).  Assume  that  the 
parameters  /3o,  /?i,  . . .,  /3P  are  precomputed. 

1.  Solve  equation  (3.19)  for  h  by  forward  substitution. 

2.  Solve  equation  (3.23)  and  form  W  via  (3.24). 

3.  Compute  R~l  by  backward  substitution,  and  form  matrix  B. 

4.  Solve  equation  (3.21)  for  g  using  a  Toeplitz  type  method. 

5.  Calculate  the  solution  vector  y  of  (3.15). 

6.  Solve  equation  (3.27)  for  r. 

7.  Form  V  via  (3.28). 

8.  Solve  (3.29)  for  s. 

9.  Compute  the  solution  vector  x  via  (3.26). 

endalgorithm 

4.  Block  Circulant  Systems 

A  block  matrix  M  =  circ(Mo,  Ai), . . . ,  Mm_i),  where  each  of  the  blocks  Mj  is  itself  an  n-th 
order  circulant,  is  called  block  circulant  with  circulant  blocks.  In  this  section  we  will  consider  the 
system  of  the  form 

(4.1)  Mx  =  6, 


which  arises  in  many  applications.  An  important  one  is  the  solution  of  elliptic  equation  with 


periodic  boundary  conditions  by  finite  difference  methods.  The  method  proposed  here  to  solve 
(4.1)  is  based  on  the  use  of  the  fast  Fourier  transform,  and  on  the  fact  that  the  coefficient  matrices 
of  the  resulting  linear  systems  are  circulants,  hence  the  methods  described  in  section  2  and  3  may 


be  employed  directly. 


$ 


,y., - .  -iVvn'r  •  >7 


•r.:”r  -t-yr : 


Denote  the  eigenvalues  of  the  n-by-n  circulant  matrices  Mj  by  . . . ,  Xn\  then  we  have 

(4.2)  Mj  =  F  DjF* 
where 

(4.3)  Dj  =  diag(X^\  . . . ,  A^),  j  =  0, 1, . . .  ,m  -  1. 

It  follows  that  the  block  circulant  matrix  M  can  be  decomposed  as 

(4.4)  M  =  QDQ\ 
where 

(4.5)  Q  =  diag(F  ,F,...,F) 
is  an  unitary  block  diagonal  matrix  of  order  (mn),  and 

D  =  circ(Do,  D\, . .  .,Dm- j) 


is  a  block  circulant. 

After  certain  exchanges  of  the  rows  and  the  columns  of  D  we  are  led  to 

(4.6)  M  =  QPdiag(N\ ,  JV2, . . . ,  Nn)PQ', 
where 

(4.7)  Nk  =  drc(xf\  X™ . X^),  k  =  1,2, . . .  ,n 

are  m-by-m  circulant  matrices  with  the  eigenvalues  of  Mj  as  their  elements,  and  P  is  an  (mn)- 
by- (mn)  permutation  matrix.  It  is  easily  seen  from  (4.6)  that  M  is  nonsingular  if  and  only  if  the 
matrices  Nk  are  nonsingular  for  all  k.  In  this  case  we  have 

(4.8)  A/-1  =  QPdiag(N^,N^\...,N~l)PQ\ 
and  the  solution  to  (4.1)  is  given  by 


(4.9) 


x  =  QPdiag(Nr\N;\ . . . ,  N^)PQ'b. 


If  at  least  one  of  the  diagonal  blocks  is  singular,  then  the  inverse  of  M  does  not  exist  and  the  least 
squares  solution  of  (4.1)  is 


(4.10) 


x  =  QPdiag(N+,N+,. . . ,  N+)PQ'b, 


which  will  be  the  solution  to  (4.1)  if  the  system  is  consistent. 

To  use  (4.9)  and  (4.10)  to  compute  the  solution  of  (4.1),  we  need  to  do  some  exchanges  of 
the  components  of  some  vectors  in  addition  to  performing  the  FFT’s  and  solving  the  circulant 
systems.  These  exchanges  can  easily  be  handled  by  partitioning  the  vector  6  and  x  into  subvectors 
to  conform  with  the  block  structure  of  the  matrix  M  as 


/  bM  \ 

( x(1)  A 

(4.11) 

6  = 

6(2) 

,  x  = 

X(2) 

V  6(m) ) 

U(m)y 

where 

(4.12) 


x( j)  =  {x\j,X2j,...,Xnj)T , 
are  subvectors  of  length  n.  Now  let 

6  (>')  =  j  =  1,2,  ...,m, 


j  =  1,2,...,  m 


(4.13) 


and  the  components  of  fcw)  are  denoted  as  in  (4.12).  Furthermore  denote 

/  ^(1)  \ 

6(2) 


(4.14) 


6  = 


\Um)J 


where 

(4.15) 

then  we  have 


=  (hi,h2,...,bkm)T,  k=  l,2,...,n, 


6  =  PQ'b. 


Let 

(4.16) 


=  ( Xkl,Xk2,...,Xkm)T 


denote  the  solution  vectors  of  the  systems 


(4.17)  Nky  =  Uk\  k  =  1,2, 
then 

(4.18)  =  Fx^, 
where 

(4.19)  x(j)  =  (xxy,  x2y , . . . ,  xny)T,  J  =  1, 2, . . . ,  m 

and  we  finaly  obtain  the  solution  of  (4.1)  by  composing  the  vector  x  using  the  subvectors  x^  via 
(4.11).  Thus  we  have  the  following  algorithm. 

Algorithm  BLCS  (BLock  Circulant  Solver)  solves  block  circulant  system  (4.1). 

1.  Transform  =  F*b^\  j  =  1,2 , . . . ,  m. 

2.  Compute  the  eigenvalues  of  My,  j  =  1,2 , . . . ,  m. 

3.  Solve  equations  (4.17)  using  algorithm  CIRS  for  all  k. 

4.  Transform  x k‘)  —  Fx^\  j  =  1,2, . . . ,  m. 

5.  Assemble  the  vectors  x^)  via  (4.11)  to  obtain  the  solution  vector  x. 

endalgorithm 

5.  Applications  to  elliptic  equations 

We  will  now  apply  the  Fourier-circulant  method  to  certain  elliptic  equations.  It  is  well  known  [5] 
that  the  approximation  of  Poisson’s  equation  on  a  rectangle  subject  to  periodic  boundary  conditions 
in  both  directions  by  the  standard  five-point  difference  scheme  on  a  uniform  mesh  results  in  the 
block  circulant  linear  system  of  the  form 

(5.1)  Mpu  =  6, 

where 

Mp  =  circ(A,  -1,0, ...,  0,-7), 

and 

A  =  circ(  4,  -1,0, ...  ,0,  -1). 


We  partition  the  unknown  vector  u  as  in  (4.11)  and  denote  its  components  as  in  (4.12).  Since 
A  is  symmetric  circulant  tridiagonal  matrix  its  eigenvalues  are,  from  (3.2), 
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A*  =  4  -  2cos(A:  -  1)  — ,  k  =  1,2,. . .  ,n. 
n 

Assume  that  n  is  even.  We  rewrite  the  eigenvalues,  by  reordering  them,  as  follows 

r  *1  =  2, 

(5.2)  |  A2  =  6, 

I  2;r  ft 

I  A2„-i  =  A2„  =  4  -  2 cos(i/ -  1)  — ,  i/ =  2,3,...,-. 
v  n  z 

The  elements  of  the  matrix  of  the  corresponding  eigenvectors(denoted  by  F  also)  are  real  and  given 

by 

=  yi 

/.\2«/-i  =  ^ cos  J(i/  -  l)(i  -  1)^  , 
fi,2u  =  sin  ^  ~  1)(*  “  1)~  >  v-2,3,...,|. 

*  =  l,2,...,n. 

In  this  situation  we  have 

(5.4)  N/c  =  circ(Afc,— 1,0, ...,  0,-1),  k  =  1,2, . . . ,  n. 

All  the  diagonal  blocks  iV*  except  Ni, which  is  singular  since  it  has  an  eigenvalue  equal  to  zero,  are 
nonsingular,  and  therefore  M  is  singular.  From  (4.10)  the  least  squares  solution  to  (5.1)  is  given 
by 

(5.5)  «  =  QPdiag(N+,Njl . N~l)PQTb , 

where  the  orthogonal  matrix  Q  is  as  in  (4.5).  The  vector  u  computed  via  (5.5)  is  the  solution  to 
(5.1)  only  if  eTb  =  0,  where  e  =  (1,1,...,  l)r,  since  e  is  the  only  null  vector  of  matrix  Mp. 

To  compute 


(5.6) 


«(»)  =  JV+fcW, 


we  use  the  algorithm  CIRS.  In  this  case  the  eigenvalues  of  Ni  are  (Assume  that  m  is  also  even.) 


'  Mi  =  0, 

M2  =  4, 

M2*/— l  =  M2t/  =  2  -  2  cos(i/  -  1)  —  ,  v  =  2, 3, . . . , 

mi 

and  the  elements  of  the  orthogonal  matrix  of  the  corresponding  eigenvectors  are 


*•* = 

<;,-,2^-i  =  y^cos  [("  -  i)(*  ~  . 

guv  =  \Z^sin  [(^  -  i)(*  “  ’ 


^  —  2, 3, . . . ,  — , 


If  we  denote  this  matrix  by 


then  (5.6)  becomes 


i  =  l,2,...,m. 


(5.7)  «(1)  =  Gdiag{0,ti2  \ . .  •  ,Mm)<Sr6(1), 

which  can  be  computed  by  real  transform  and  its  inverse. 

The  other  equations 

(5.8)  Nk  «<*>  =  ¥k) 

are  nonsingular  circulant  tridiagonal,  and  the  coefficient  matrices  have  the  forms 

/At  -1  -1\ 


-1  A*  -1 


,  k  =  2,3,. ..,n. 


-1  A*  -1 
-1  A*. 


Since  the  eigenvalues  satisfy  the  relation  At  >  2  for  k  —  2, 3, . . . ,  n,  the  quadratic  equations 

(5.9)  r2  -  Atr  +  1  =  0 


f.-TV'A; 


k  =  2,3 , . . . ,  n, 


and  ei  and  em  are  the  first  and  the  last  column  of  the  m-th  order  identity  matrix.  Thus  the 
algorithm  proceeds  as  follows. 

Algorithm  PPS  (Periodic  Poisson  Solver)  solves  block  circulant  system  (5.1). 

1.  Transform  6^')  =  F*b j  =  1,2, . . . ,  m. 

2.  Compute  the  eigenvalues  of  A  via  (5.2). 

3.  Calculate  via  (5.7). 

4.  Solve  equation  (5.8)  for  k  =  2, 3, . . . ,  n  using  algorithm  CTS. 

5.  Transform  x (*’)  =  Fx^\  j  =  1,2, . . . ,  m. 

6.  Assemble  the  vectors  x W  via  (4.11)  to  obtain  the  solution  vector  x. 

endalgorithm 

The  complexity  of  the  algorithm  is  0(mn  log2  n).  Although  it  is  the  same  as  cyclic  reduction, 
the  algorithm  is  more  efficient  than  cyclic  reduction  since  it  is  a  FFT-based  method(see  [16]) 

Another  kind  of  elliptic  equations  to  which  the  method  may  be  applied  is  the  biharmonic 
equation 

(5.13)  V4u(x,y)  =  (x,y)€fi 

subject  to  periodic  conditions. 

When  the  region  Q  is  a  rectangle,  and  the  mesh  lengths  in  both  directions  are  the  same,  and 
the  standard  13-point  difference  approximation  of  the  biharmonic  operator  is  used,  the  resulting 
system  of  linear  equations  is  of  the  form 

(5.14)  Mtu  =  5, 


where 


and 


iWi  =  circ(C,  B,J,0,...,0, 1,  B), 

C  —  circ( 20,  -8, 1,0, ..  .,0, 1,  -8), 
B  =  c»rc(-8,2,0,  ...,0,2). 


Since  C  and  B  are  symmetric  their  eigenvalues  are  real  and  given  by 
Ai=6, 


(5.15) 


A2  =  38, 

A21/-1  =  A2„  =  20  -  16cos(y  —  1)  —  +  2cos2(i/  —  1)  —  ,  v 

n  n 
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and 

(5.16) 

respectively.  Thus 


f  m  =  -4, 

P2  =  -12, 

„  „  ,  „ .  2ir  „  „  m 

P2v-i  =  P2i>  =  -8  +  4cos(i/-  1)  —  ,  u  =  2,3,...,—, 

771  2 


(5.17)  JVt  =  «>c(A*,m*,1,0,...,0,1,/i*),  k  =  1,2, ...  ,n. 

The  system 

(5.18)  =  6*1) 

is  singular,  and  can  be  solved  by  using  the  algorithm  CIRS,  and  the  others 

(5.19)  jVfctiW  =  g(*>,  fc  =  2,3,..  .,ti 

are  circulant  pentadiagonal,  and  can  be  solved  by  the  method  of  §3.  In  this  case,  p  —  2  and  the 
factor  1(2),  which  is  quadratic,  may  be  found  by  solving  the  equations(  see  [7]) 

(5.20)  p 4  +  pkp3  +  \kp2  +  pkp  +  1  =  0. 

Since  p\  —  4A*  +  8  =  0  for  k  =  2,3, . . . ,  71,  the  equations  (5.20)  have  two  real  roots  that  each 
has  multiplicity  2,  and  they  are  given  by 


(5.21) 
and 

(5.22) 

It  is  easily  seen  that 

and 


Pi 


(*)  _  i 

1 


-pk  +  \Jp\  -  16 


P2 


(*)  =  I 


16 


p[k)  >  1, 


p[k)  <  1, 


0W  -  — 

Pl  ~~  n(k) 

P  2 

since  pk  <  —4  for  k  —  2, 3, . . .  ,71,  and  the  factor  lk{z)  for  the  matrix  Nk  is  found  to  be 

lk(z)=0{ok)+0[k)z  +  0{2k)z2, 

where 

=  -2, 


(5.23) 


gW  =  — 

-  (*)  ■ 
Pl 
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Thus  we  have  the  following  0(mn  log2  n)  algorithm. 

Algorithm  PBS  (Periodic  Biharmonic  Solver)  solves  block  circulant  system  (5.14). 

1.  Transform  6^’)  =  F*b(]\  j  =  1, 2, . . . ,  m. 

2.  Compute  the  eigenvalues  of  A  and  B  via  (5.15)  and  (5.16),  respectively. 

3.  Solve  equation  (5.18)  for  using  algorithm  CIRS. 

4.  For  k  =  2, 3, . . . ,  n, 

4.1.  calculate  /3^k\  fi[k\  via  (5.21)  and  (5.23),  and 

4.2.  solve  equations  (5.19)  for  using  algorithm  BCS. 

5.  Transform  =  Fx^\  j  =  1, 2, . . . ,  m. 

6.  assemble  the  vectors  via  (4.11)  to  obtain  the  solution  vector  x. 

endalgorithm 


6.  Numerical  experiments 

The  methods  described  in  the  previous  section  were  tried  on  the  FPS-164  of  the  Department  of 
Computer  Science,  Yale  University.  The  programs  were  written  and  timed  in  FORTRAN,  and  the 
FFT  subroutines  used  are  part  of  the  NA  Packege  of  the  Research  Center  for  Scientific  Computation, 
Yale  University. 

In  each  test  case,  a  periodic  function  was  arbitrarily  chosen  and  its  values  at  the  mesh  points  of 
an  N  X  N  grid,  which  are  considered  to  be  the  “true”  solution  of  (5.1)  and  (5.14),  were  calculated, 
and  the  corresponding  right  hand  sides  were  computed  with  the  use  of  (5.1)  and  (5.14),  respectively. 
Approximate  solutions  to  (5.1)  and  (5.14)  were  then  generated  by  using  the  algorithms  described 
in  §5  as  well  as  the  cyclic  reduction  for  Poisson’s  equation.  The  maximun  errors  of  all  tests  are  less 
than  10-10,  indicating  that  the  algorithms  presented  in  this  paper  are  stable.  The  execution  time 
for  each  case  are  summarized  in  Table  1.  It  should  be  noted  that  the  Foureir-circulant  method  is 
more  efficient  than  cyclic  reduction. 


N 

Algorithm  PPS 

Cyclic  reduction 

Algorithm  PBS 

16 

0.246  x  10- 1 

0.816  x  10_1 

0.287  x  10- 1 

32 

0.696  x  10" 1 

|—  ||  |  \j  || 

0.815  x  10-i 

64 

0.258  x  10° 

0.155  X  10l 

128 

■HQQH 

0.105  X  101 

Table  1:  Execution  times  (millisecond)  on  FPS-164 
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